% *** MRSI 2D Reconstruction *** %
%   Author : Sungjin Kim
%   Created Date : 23 Oct, 2013
%   Modified Date : 23 Oct, 2013
%   Version : 0.0
%   Description : Reconstruct 2D MRSI from single channel FID dataset.
%                   
% *********************************************************************** %

function [mrsi, headerinfo] = mrsi2Drecon(filename, cell_args)
    % get header information.
    cell_args = {3}; % default value; display nothing
    [~, headerinfo] = readfidn(filename, 1, 1, cell_args);

    nblocks=headerinfo{1, 1}(1) % header info
    ntraces=headerinfo{1, 1}(2) % status info
    dataP=headerinfo{1, 1}(3)/2; % data points
    
    % get data from blocks and traces
    complexFID=cell(nblocks, ntraces);

    for ix_trace=1:1:ntraces
        for ix_block=1:1:nblocks
            [complexFID{ix_block, ix_trace}, ~] = readfidn(filename, ix_block, ix_trace, cell_args);
        end
    end

    % 1st FFT - to convert fid into spectra
    cmplxSpec = cell(nblocks, ntraces);
    for ix_trace=1:1:ntraces
        for ix_block=1:1:nblocks
            cmplxSpec{ix_block, ix_trace} = fft(complexFID{ix_block, ix_trace});
            cmplxSpec{ix_block, ix_trace} = circshift(cmplxSpec{ix_block, ix_trace}, dataP/2);
        end
    end
    
    % convert cell dataset into 3D array
    arr_cmplx2Dmrsi = zeros(nblocks, ntraces, dataP);
    rot_arr_cmplx2Dmrsi = zeros(ntraces, nblocks, dataP);
    
    for indx_block = 1:1:nblocks 
        for indx_trace = 1:1:ntraces
            for indx_data = 1:1:dataP
                arr_cmplx2Dmrsi(indx_block, indx_trace, indx_data) = cmplxSpec{indx_block, indx_trace}(indx_data);                
            end
        end
    end
    
    for indx_data = 1:1:dataP
        rot_arr_cmplx2Dmrsi(:, :, indx_data) = rot90(arr_cmplx2Dmrsi(:, :, indx_data), -1);
    end

    % display k-space
    figure(11);
    im=sum(abs(rot_arr_cmplx2Dmrsi),3);
    imagesc(im);
    title('K-Space');
    colormap('gray');
        
    % 2D fft
    mrsi = fft2(rot_arr_cmplx2Dmrsi);
    mrsi = circshift(mrsi, [ntraces/2, nblocks/2]);
    
    return;
end % end of function